To Do:

  • Fill missing data
  • Update Pinery Climate (96 - present; Thedford & our data)
  • Growing season T (May - Sept)
  • Vapour pressure defecit (growing season T & previous Aug-Oct, Williams et al.)

1 Packages

library(ggplot2) #for plotting
library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ tibble  3.0.1     ✓ dplyr   1.0.0
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ✓ purrr   0.3.4
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(dplyr)
library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

2 Read in Data

The Pinery dataset was compiled in excel from monthly mean temperature and total precipitation data using open access Environment Canada records for stations within 30km of the Pinery with at least 10 consecutive years of data (http://climate.weather.gc.ca/historical_data/search_historic_data_e.html).

The London datasets are homogenized data from Environment Canada that were converted to tidy format in excel.

pinery <- read.csv("data/pinery_climate.csv")
pinery$station <- as.factor(pinery$station) #Make a factor 
str(pinery)
## 'data.frame':    2931 obs. of  5 variables:
##  $ year        : int  1882 1882 1882 1882 1882 1882 1883 1883 1883 1883 ...
##  $ month       : int  7 8 9 10 11 12 4 5 6 7 ...
##  $ mean_temp   : num  NA NA 16.1 11.4 2.9 -3.9 NA NA NA NA ...
##  $ total_precip: num  15.5 182.6 57.9 41.7 28.7 ...
##  $ station     : Factor w/ 9 levels "Arkona","Centralia A",..: 1 1 1 1 1 1 8 8 8 8 ...
london_precip <- read.csv("data/london_precip.csv")
london_temp <- read.csv("data/london_temp.csv")
london_precip$month <- as.factor(london_precip$month) #Make a factor 
london_temp$month <- as.factor(london_temp$month) #Make a factor 
climateNA <- read_csv("data/climateNA_monthly.csv") 
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_double(),
##   Year = col_double(),
##   Month = col_double(),
##   Tmax = col_double(),
##   Tmin = col_double(),
##   Tave = col_double(),
##   PPT = col_double(),
##   Rad = col_double(),
##   Eref = col_double(),
##   CMD = col_double(),
##   RH = col_double()
## )
climateNA.temp <- climateNA %>% select(Year, Month, Tave)
climateNA.precip <- climateNA %>% select(Year, Month, PPT)

3 Create Regional Data

3.1 Clean Data

Here I go through each station for the temp and precip data, removing years where all 12 months aren’t present so that I can later calculate the mean values for each year. I did this manually a while ago and there is probably a more efficient way to do it, however it works and it’s already done so I left it as is.

#Arkona
Arkona <- pinery %>% filter(station == "Arkona") #Filter out station of interest

Arkona_T <- Arkona %>% select(-total_precip) %>% filter(!is.na(Arkona$mean_temp)==T) #remove precipitation data and "NA"
table(Arkona_T$year) #Generate table after cleaning
## 
## 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895 1896 1897 
##    4   11   11   12   12   12   12   12   12   12   12    5   12   12   12   12 
## 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 
##   12   12   10   12    4   12   12   12   12   12   12   12   12   12   12   12 
## 1914 1915 
##   12    3
Arkona_T2 <- Arkona_T %>% filter(year != 1882 & year != 1883 & 
                                   year != 1884 & year != 1893 & year != 1900 & year != 1902 & year != 1915) #remove years with < 12 months
table(Arkona_T2$year) #Generate table after cleaning
## 
## 1885 1886 1887 1888 1889 1890 1891 1892 1894 1895 1896 1897 1898 1899 1901 1903 
##   12   12   12   12   12   12   12   12   12   12   12   12   12   12   12   12 
## 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914 
##   12   12   12   12   12   12   12   12   12   12   12
Arkona_P <- Arkona %>% select(-mean_temp) %>% filter(!is.na(Arkona$total_precip)==T) #remove temperature data and "NA"
#table(Arkona_P$year) #show how many observations in each year
Arkona_P2 <- Arkona_P %>% filter(year != 1882 & year != 1893 & year != 1902 & year != 1915)
#table(Arkona_P2$year) #Generate table after cleaning

#Centralia A
CentraliaA <- pinery %>% filter(station == "Centralia A")

CentraliaA_T <- CentraliaA %>% select(-total_precip) %>% filter(!is.na(CentraliaA$mean_temp)==T) #remove "NA"
#table(CentraliaA_T$year)  
CentraliaA_T2 <- CentraliaA_T %>% filter(year != 1945 & year != 1946 & year != 1947 & year != 1966)
#table(CentraliaA_T2$year)  

CentraliaA_P <- CentraliaA %>% select(-mean_temp) %>% filter(!is.na(CentraliaA$total_precip)==T) #remove "NA"
#table(CentraliaA_P$year)  
CentraliaA_P2 <- CentraliaA_P %>% filter(year != 1946 & year != 1947 & year != 1966)
#table(CentraliaA_P2$year)  

#Exeter
Exeter <- pinery %>% filter(station == "Exeter")

Exeter_T <- Exeter %>% select(-total_precip) %>% filter(!is.na(Exeter$mean_temp)==T) #remove "NA"
#table(Exeter_T$year)  
Exeter_T2 <- Exeter_T %>% filter(year != 1967 & year != 1976 & year != 1985)
#table(Exeter_T2$year)  

Exeter_P <- Exeter %>% select(-mean_temp) %>% filter(!is.na(Exeter$total_precip)==T) #remove "NA"
#table(Exeter_P$year)  
Exeter_P2 <- Exeter_P %>% filter(year != 1961 & year != 1962 & year != 1964 & year != 1976 & year != 1985)
#table(Exeter_P2$year) 

#Forest
Forest <- pinery %>% filter(station == "Forest")

Forest_T <- Forest %>% select(-total_precip) %>% filter(!is.na(Forest$mean_temp)==T) #remove "NA"
#table(Forest_T$year)  
Forest_T2 <- Forest_T %>% filter(year != 1924 & year != 1940 & year != 1962 & year != 1963)
#table(Forest_T2$year) 

Forest_P <- Forest %>% select(-mean_temp) %>% filter(!is.na(Forest$total_precip)==T) #remove "NA"
#table(Forest_P$year)  
Forest_P2 <- Forest_P %>% filter(year > 1933 & year != 1937 & year!= 1943 & year != 1955 & year != 1960 & year < 1960)
#table(Forest_P2$year) 

#Lucan
Lucan <- pinery %>% filter(station == "Lucan")

Lucan_T <- Lucan %>% select(-total_precip) %>% filter(!is.na(Lucan$mean_temp)==T) #remove "NA"
#table(Lucan_T$year)  
Lucan_T2 <- Lucan_T %>% filter(year != 1915 & year != 1924 & year != 1948 & year != 1949 & year != 1950
                               & year != 1951 & year != 1952 & year != 1953 & year != 1954 & year != 1955
                               & year != 1956 & year != 1957 & year != 1958 & year != 1959 & year != 1961)
#table(Lucan_T2$year) 

Lucan_P <- Lucan %>% select(-mean_temp) %>% filter(!is.na(Lucan$total_precip)==T) #remove "NA"
#table(Lucan_P$year)  
Lucan_P2 <- Lucan_P %>% filter(year != 1915 & year != 1923 & year != 1943 & year < 1948)
#table(Lucan_P2$year)

#Pinery
Pinery <- pinery %>% filter(station == "Pinery")

Pinery_T <- Pinery %>% select(-total_precip) %>% filter(!is.na(Pinery$mean_temp)==T) #remove "NA"
#table(Pinery_T$year)  
Pinery_T2 <- Pinery_T %>% filter(year != 1979 & year != 1984)
#table(Pinery_T2$year) 

Pinery_P <- Pinery %>% select(-mean_temp) %>% filter(!is.na(Pinery$total_precip)==T) #remove "NA"
#table(Pinery_P$year)  
Pinery_P2 <- Pinery_P %>% filter(year != 1979 & year != 1984)
#table(Pinery_P2$year)

#Strathrow-Mullifary
Strathroy <- pinery %>% filter(station == "Strathroy-Mullifary")

Strathroy_T <- Strathroy %>% select(-total_precip) %>% filter(!is.na(Strathroy$mean_temp)==T) #remove "NA"
#table(Strathroy_T$year)  
Strathroy_T2 <- Strathroy_T %>% filter(year != 2009 & year != 2010 & year != 2013 & year != 2014 & year != 2015 & year < 2017)
#table(Strathroy_T2$year) 

Strathroy_P <- Strathroy %>% select(-mean_temp) %>% filter(!is.na(Strathroy$total_precip)==T) #remove "NA"
#table(Strathroy_P$year)  
Strathroy_P2 <- Strathroy_P %>% filter(year != 2009 & year != 2010 & year != 2012 & year != 2013 
                                       & year != 2014 & year != 2015 & year < 2017)
#table(Strathroy_P2$year)

#Thedford (no temp)
Thedford <- pinery %>% filter(station == "Thedford")

Thedford_P <- Thedford %>% select(-mean_temp) %>% filter(!is.na(Thedford$total_precip)==T) #remove "NA"
#table(Thedford_P$year)  
Thedford_P2 <- Thedford_P %>% filter(year != 1883 & year != 1897)
#table(Thedford_P2$year)

#Watford (can't use precipitation data)
Watford <- pinery %>% filter(station == "Watford")

Watford_T <- Watford %>% select(-total_precip) %>% filter(!is.na(Watford$mean_temp)==T) #remove "NA"
#table(Watford_T$year)  
Watford_T2 <- Watford_T %>% filter(year == 1927 | year == 1928)
#table(Watford_T2$year) 

Watford_P <- Watford %>% select(-mean_temp) %>% filter(!is.na(Watford$total_precip)==T) #remove "NA"
#table(Watford_P$year)  

Re-combine Data

pinery_temp <- rbind(Arkona_T2, CentraliaA_T2, Exeter_T2, Forest_T2, Lucan_T2, Pinery_T2, Strathroy_T2, Watford_T2)
pinery_precip <- rbind(Arkona_P2, CentraliaA_P2, Exeter_P2, Forest_P2, Lucan_P2, Pinery_P2, Strathroy_P2, Thedford_P2)

Ensure the cleaning worked by plotting histograms of the number of observations (months) in each year (should only be multiples of 12)

#histogram (temperature)
graph <- ggplot(pinery_temp, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
  geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
    geom_hline(yintercept = 24, linetype="solid", color = "red") + 
    geom_hline(yintercept = 36, linetype="solid", color = "red") + 
  theme_bw()
graph

#histogram (precipitation)
graph <- ggplot(pinery_precip, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
   geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
    geom_hline(yintercept = 24, linetype="solid", color = "red") + 
  theme_bw()
graph

3.2 Calculate Average (month-year)

Calculate one mean temperature and total precipitation for each year/month where more than one station of data is available

#Temperature
str(pinery_temp)
## 'data.frame':    1944 obs. of  4 variables:
##  $ year     : int  1885 1885 1885 1885 1885 1885 1885 1885 1885 1885 ...
##  $ month    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ mean_temp: num  -9.4 -13.5 -8.9 4.2 12.1 16.4 21.1 17 14.8 8 ...
##  $ station  : Factor w/ 9 levels "Arkona","Centralia A",..: 1 1 1 1 1 1 1 1 1 1 ...
pinery_temp_mean <- pinery_temp %>%
  group_by(year,month) %>%
  summarise(mean_temp_year = mean(mean_temp), sd_temp = sd(mean_temp), n = n()) %>%
  mutate(se_temp = sd_temp / sqrt(n), lower = (mean_temp_year - se_temp), upper = (mean_temp_year + se_temp)) %>% round(2)
## `summarise()` regrouping output by 'year' (override with `.groups` argument)
pinery_temp_mean$month <- as.factor(pinery_temp_mean$month)

#Precipitation
pinery_precip_mean <- pinery_precip %>%
  group_by(year,month) %>%
  summarise(total_precip_year = mean(total_precip), sd_precip = sd(total_precip), n = n()) %>%
  mutate(se_precip = sd_precip / sqrt(n), lower = total_precip_year - sd_precip, upper = total_precip_year + sd_precip) %>% round(2)
## `summarise()` regrouping output by 'year' (override with `.groups` argument)
head(pinery_precip_mean)
## # A tibble: 6 x 8
## # Groups:   year [1]
##    year month total_precip_year sd_precip     n se_precip lower upper
##   <dbl> <dbl>             <dbl>     <dbl> <dbl>     <dbl> <dbl> <dbl>
## 1  1883     1              30.2        NA     1        NA    NA    NA
## 2  1883     2              94.7        NA     1        NA    NA    NA
## 3  1883     3              70.6        NA     1        NA    NA    NA
## 4  1883     4              33          NA     1        NA    NA    NA
## 5  1883     5             112.         NA     1        NA    NA    NA
## 6  1883     6             158.         NA     1        NA    NA    NA
ggplot(pinery_precip_mean, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey")

ggplot(pinery_temp_mean, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey")

3.3 Fill Empty Years

3.3.1 Temperature

Fill empty months-years with data from London

p.t <- pinery_temp_mean %>% select(year, month, mean_temp_year)

l.t <- london_temp %>% filter(month == "Jan" |
                                month == "Feb" |
                                month == "Mar" |
                                month == "Apr" |
                                month == "May" |
                                month == "Jun" |
                                month == "Jul" |
                                month == "Aug" |
                                month == "Sep" |
                                month == "Oct" |
                                month == "Nov" |
                                month == "Dec") %>% 
  mutate(month = recode(month, "Jan"="1", 
                                     "Feb"="2", 
                                     "Mar"="3",
                                     "Apr"="4",
                                     "May"="5", 
                                     "Jun"="6", 
                                     "Jul"="7",
                                     "Aug"="8",
                                     "Sep"="9",
                                     "Oct"="10",
                                     "Nov"="11",
                                     "Dec"="12"))
PL.temp <- full_join(l.t, p.t) %>%
  rename("pinery" = "mean_temp_year") %>%
  rename("london" = "mean_temp")
## Joining, by = c("year", "month")
head(PL.temp)
##   year month london pinery
## 1 1883     1     NA     NA
## 2 1884     1  -10.1     NA
## 3 1885     1   -9.2   -9.4
## 4 1886     1   -7.2   -7.2
## 5 1887     1   -8.3   -8.4
## 6 1888     1  -10.1   -9.6

Select all rows with regional data - source = Pinery Select all rows with NA - source = London

PL.1 <- PL.temp %>%  #use london value where pinery is absent
  filter(is.na(pinery)) %>%
  mutate(temp = london)
PL.1$source <- "London"

PL.2 <- PL.temp %>% # use pinery value where pinery is present
  filter(!is.na(pinery)) %>%
  mutate(temp = pinery)
PL.2$source <- "Pinery"

PL.temp2 <- rbind(PL.1, PL.2) %>%
  arrange(year)

(nrow(PL.2) + nrow(PL.1)) == nrow(PL.temp) #check equal observations
## [1] TRUE
PL.temp2 %>% filter(is.na(temp))
##   year month london pinery temp source
## 1 1883     1     NA     NA   NA London
## 2 1883     2     NA     NA   NA London
## 3 2014     9     NA     NA   NA London

Fill leftover NAs with climateNA values

head(climateNA.temp)
## # A tibble: 6 x 3
##    Year Month  Tave
##   <dbl> <dbl> <dbl>
## 1  1901     1  -4.2
## 2  1902     1  -5.6
## 3  1903     1  -5.1
## 4  1904     1  -9.5
## 5  1905     1  -7.5
## 6  1906     1  -0.7
climateNA.temp <- climateNA.temp %>%
  rename("month" = "Month") %>%
  rename("year" = "Year") %>%
  rename("climateNA" = "Tave")

climateNA.temp$month <- as.factor(climateNA.temp$month)

PL.temp3 <- full_join(PL.temp2, climateNA.temp)
## Joining, by = c("year", "month")
PL.temp3 <- PL.temp3[, c(1:4,7,5,6)] #re-order columns

head(PL.temp3)
##   year month london pinery climateNA temp source
## 1 1883     1     NA     NA        NA   NA London
## 2 1883     2     NA     NA        NA   NA London
## 3 1883     3   -6.5     NA        NA -6.5 London
## 4 1883     4    5.1     NA        NA  5.1 London
## 5 1883     5   10.0     NA        NA 10.0 London
## 6 1883     6   17.3     NA        NA 17.3 London
PL.3 <- PL.temp3 %>%  
  filter(is.na(temp)) %>%
  mutate(temp = climateNA)
PL.3$source <- "ClimateNA"

PL.4 <- PL.temp3 %>% 
  filter(!is.na(temp))

PL.temp4 <- rbind(PL.3, PL.4) %>%
  arrange(year) %>% 
  filter(year > 1890) #missing values from 1890s remain

(nrow(PL.3) + nrow(PL.4)) == nrow(PL.temp3) #check equal observations
## [1] TRUE
PL.temp4 %>% filter(is.na(temp)) #values still missing
## [1] year      month     london    pinery    climateNA temp      source   
## <0 rows> (or 0-length row.names)
ggplot(PL.temp4, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
  geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
  theme_bw()

3.3.2 Precipitation

p.p <- pinery_precip_mean %>% select(year, month, total_precip_year)

l.p <- london_precip %>% filter(month == "Jan" |
                                month == "Feb" |
                                month == "Mar" |
                                month == "Apr" |
                                month == "May" |
                                month == "Jun" |
                                month == "Jul" |
                                month == "Aug" |
                                month == "Sep" |
                                month == "Oct" |
                                month == "Nov" |
                                month == "Dec") %>% 
  mutate(month = recode(month, "Jan"="1", 
                                     "Feb"="2", 
                                     "Mar"="3",
                                     "Apr"="4",
                                     "May"="5", 
                                     "Jun"="6", 
                                     "Jul"="7",
                                     "Aug"="8",
                                     "Sep"="9",
                                     "Oct"="10",
                                     "Nov"="11",
                                     "Dec"="12"))
p.p$month <- as.factor(p.p$month)
PL.precip <- full_join(l.p, p.p) %>%
  rename("pinery" = "total_precip_year") %>%
  rename("london" = "total_precip")
## Joining, by = c("year", "month")
PL.5 <- PL.precip %>%  #use london value where pinery is absent
  filter(is.na(pinery)) %>%
  mutate(precip = london)
PL.5$source <- "London"

PL.6 <- PL.precip %>% # use pinery value where pinery is present
  filter(!is.na(pinery)) %>% 
  mutate(precip = pinery)
PL.6$source <- "Pinery"

PL.precip2 <- rbind(PL.5, PL.6) %>%
  arrange(year)

(nrow(PL.5) + nrow(PL.6)) == nrow(PL.precip2) #check equal observations
## [1] TRUE
PL.precip2 %>% filter(is.na(precip)) #values still missing
##    year month london pinery precip source
## 1  2009     4     NA     NA     NA London
## 2  2009     5     NA     NA     NA London
## 3  2009     6     NA     NA     NA London
## 4  2009     7     NA     NA     NA London
## 5  2009     8     NA     NA     NA London
## 6  2009     9     NA     NA     NA London
## 7  2009    10     NA     NA     NA London
## 8  2010     4     NA     NA     NA London
## 9  2010     5     NA     NA     NA London
## 10 2010     6     NA     NA     NA London
## 11 2010     7     NA     NA     NA London
## 12 2010     8     NA     NA     NA London
## 13 2010     9     NA     NA     NA London
## 14 2010    10     NA     NA     NA London
## 15 2012     4     NA     NA     NA London
## 16 2012     5     NA     NA     NA London
## 17 2012     6     NA     NA     NA London
## 18 2012     7     NA     NA     NA London
## 19 2012     8     NA     NA     NA London
## 20 2012     9     NA     NA     NA London
## 21 2012    10     NA     NA     NA London
## 22 2013     4     NA     NA     NA London
## 23 2013     5     NA     NA     NA London
## 24 2013     6     NA     NA     NA London
## 25 2013     7     NA     NA     NA London
## 26 2013     8     NA     NA     NA London
## 27 2013     9     NA     NA     NA London
## 28 2013    10     NA     NA     NA London
## 29 2014     4     NA     NA     NA London
## 30 2014     5     NA     NA     NA London
## 31 2014     6     NA     NA     NA London
## 32 2014     7     NA     NA     NA London
## 33 2014     8     NA     NA     NA London
## 34 2014     9     NA     NA     NA London
## 35 2014    10     NA     NA     NA London
## 36 2015     4     NA     NA     NA London
## 37 2015     5     NA     NA     NA London
## 38 2015     6     NA     NA     NA London
## 39 2015     7     NA     NA     NA London
## 40 2015     8     NA     NA     NA London
## 41 2015     9     NA     NA     NA London
## 42 2015    10     NA     NA     NA London
## 43 2017     4     NA     NA     NA London
## 44 2017     5     NA     NA     NA London
## 45 2017     6     NA     NA     NA London
## 46 2017     7     NA     NA     NA London
## 47 2017     8     NA     NA     NA London
## 48 2017     9     NA     NA     NA London
## 49 2017    10     NA     NA     NA London
## 50 2017    11     NA     NA     NA London
## 51 2017    12     NA     NA     NA London

Fill leftover NAs with climateNA values

head(climateNA.temp)
## # A tibble: 6 x 3
##    year month climateNA
##   <dbl> <fct>     <dbl>
## 1  1901 1          -4.2
## 2  1902 1          -5.6
## 3  1903 1          -5.1
## 4  1904 1          -9.5
## 5  1905 1          -7.5
## 6  1906 1          -0.7
climateNA.precip <- climateNA.precip %>%
  rename("month" = "Month") %>%
  rename("year" = "Year") %>%
  rename("climateNA" = "PPT")

climateNA.precip$month <- as.factor(climateNA.precip$month)

PL.precip3 <- full_join(PL.precip2, climateNA.precip)
## Joining, by = c("year", "month")
PL.precip3 <- PL.precip3[, c(1:4,7,5,6)] #re-order columns

head(PL.temp3)
##   year month london pinery climateNA temp source
## 1 1883     1     NA     NA        NA   NA London
## 2 1883     2     NA     NA        NA   NA London
## 3 1883     3   -6.5     NA        NA -6.5 London
## 4 1883     4    5.1     NA        NA  5.1 London
## 5 1883     5   10.0     NA        NA 10.0 London
## 6 1883     6   17.3     NA        NA 17.3 London
PL.7 <- PL.precip3 %>%  
  filter(is.na(precip)) %>%
  mutate(precip = climateNA)
PL.7$source <- "ClimateNA"

PL.8 <- PL.precip3 %>% 
  filter(!is.na(precip))

PL.7
##    year month london pinery climateNA precip    source
## 1  2009     4     NA     NA       132    132 ClimateNA
## 2  2009     5     NA     NA        67     67 ClimateNA
## 3  2009     6     NA     NA       108    108 ClimateNA
## 4  2009     7     NA     NA        72     72 ClimateNA
## 5  2009     8     NA     NA        90     90 ClimateNA
## 6  2009     9     NA     NA        40     40 ClimateNA
## 7  2009    10     NA     NA       116    116 ClimateNA
## 8  2010     4     NA     NA        64     64 ClimateNA
## 9  2010     5     NA     NA       103    103 ClimateNA
## 10 2010     6     NA     NA       132    132 ClimateNA
## 11 2010     7     NA     NA       100    100 ClimateNA
## 12 2010     8     NA     NA        36     36 ClimateNA
## 13 2010     9     NA     NA        76     76 ClimateNA
## 14 2010    10     NA     NA        57     57 ClimateNA
## 15 2012     4     NA     NA        58     58 ClimateNA
## 16 2012     5     NA     NA        50     50 ClimateNA
## 17 2012     6     NA     NA        52     52 ClimateNA
## 18 2012     7     NA     NA        67     67 ClimateNA
## 19 2012     8     NA     NA        70     70 ClimateNA
## 20 2012     9     NA     NA        73     73 ClimateNA
## 21 2012    10     NA     NA       128    128 ClimateNA
## 22 2013     4     NA     NA       115    115 ClimateNA
## 23 2013     5     NA     NA        71     71 ClimateNA
## 24 2013     6     NA     NA       131    131 ClimateNA
## 25 2013     7     NA     NA        91     91 ClimateNA
## 26 2013     8     NA     NA        85     85 ClimateNA
## 27 2013     9     NA     NA        51     51 ClimateNA
## 28 2013    10     NA     NA       123    123 ClimateNA
## 29 2014     4     NA     NA        81     81 ClimateNA
## 30 2014     5     NA     NA        81     81 ClimateNA
## 31 2014     6     NA     NA        88     88 ClimateNA
## 32 2014     7     NA     NA        90     90 ClimateNA
## 33 2014     8     NA     NA       103    103 ClimateNA
## 34 2014     9     NA     NA        95     95 ClimateNA
## 35 2014    10     NA     NA        84     84 ClimateNA
## 36 2015     4     NA     NA        74     74 ClimateNA
## 37 2015     5     NA     NA        95     95 ClimateNA
## 38 2015     6     NA     NA       118    118 ClimateNA
## 39 2015     7     NA     NA        65     65 ClimateNA
## 40 2015     8     NA     NA        79     79 ClimateNA
## 41 2015     9     NA     NA        59     59 ClimateNA
## 42 2015    10     NA     NA        68     68 ClimateNA
## 43 2017     4     NA     NA       118    118 ClimateNA
## 44 2017     5     NA     NA       122    122 ClimateNA
## 45 2017     6     NA     NA       116    116 ClimateNA
## 46 2017     7     NA     NA        75     75 ClimateNA
## 47 2017     8     NA     NA        82     82 ClimateNA
## 48 2017     9     NA     NA        86     86 ClimateNA
## 49 2017    10     NA     NA       107    107 ClimateNA
## 50 2017    11     NA     NA       110    110 ClimateNA
## 51 2017    12     NA     NA        50     50 ClimateNA
## 52 2018     1     NA     NA        66     66 ClimateNA
## 53 2018     2     NA     NA       107    107 ClimateNA
## 54 2018     3     NA     NA        68     68 ClimateNA
## 55 2018     4     NA     NA       102    102 ClimateNA
## 56 2018     5     NA     NA        67     67 ClimateNA
## 57 2018     6     NA     NA        97     97 ClimateNA
## 58 2018     7     NA     NA        89     89 ClimateNA
## 59 2018     8     NA     NA       180    180 ClimateNA
## 60 2018     9     NA     NA        71     71 ClimateNA
## 61 2018    10     NA     NA       100    100 ClimateNA
## 62 2018    11     NA     NA       136    136 ClimateNA
## 63 2018    12     NA     NA        58     58 ClimateNA
PL.precip4 <- rbind(PL.7, PL.8) %>%
  arrange(year) %>%
  filter(year >=1890)

(nrow(PL.7) + nrow(PL.8)) == nrow(PL.precip3) #check equal observations
## [1] TRUE
PL.precip4 %>% filter(is.na(precip)) #values still missing
## [1] year      month     london    pinery    climateNA precip    source   
## <0 rows> (or 0-length row.names)
ggplot(PL.precip4, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
  geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
  theme_bw()

4 Monthly Data

4.1 Climograph

study_t <- pinery_temp_mean %>% 
  group_by(month) %>% 
  summarize(mean_temp = mean(mean_temp_year) %>% round(2)) 
## `summarise()` ungrouping output (override with `.groups` argument)
study_t
## # A tibble: 12 x 2
##    month mean_temp
##    <fct>     <dbl>
##  1 1         -5.64
##  2 2         -5.78
##  3 3         -0.64
##  4 4          6.37
##  5 5         12.8 
##  6 6         18.2 
##  7 7         20.8 
##  8 8         19.9 
##  9 9         16.2 
## 10 10         9.83
## 11 11         3.28
## 12 12        -2.8
study_p <- pinery_precip_mean %>%
  group_by(month) %>%
  summarize(total_precip= mean(total_precip_year), sd_precip = sd(total_precip_year), n = n()) %>%
  mutate(se_precip = sd_precip / sqrt(n), lower = total_precip - se_precip, upper = total_precip + se_precip %>% round(2)) 
## `summarise()` ungrouping output (override with `.groups` argument)
study_p
## # A tibble: 12 x 7
##    month total_precip sd_precip     n se_precip lower upper
##    <dbl>        <dbl>     <dbl> <int>     <dbl> <dbl> <dbl>
##  1     1         80.5      32.6   122      2.95  77.6  83.5
##  2     2         65.4      26.8   122      2.42  63.0  67.8
##  3     3         64.4      26.0   122      2.35  62.1  66.8
##  4     4         72.5      32.0   122      2.90  69.7  75.4
##  5     5         82.0      38.0   122      3.44  78.6  85.5
##  6     6         78.5      36.3   122      3.28  75.2  81.8
##  7     7         81.4      48.6   122      4.40  77.0  85.8
##  8     8         79.3      43.4   122      3.93  75.4  83.3
##  9     9         86.2      51.0   122      4.61  81.6  90.8
## 10    10         82.9      41.0   122      3.72  79.2  86.6
## 11    11         89.1      32.2   122      2.91  86.2  92.0
## 12    12         87.2      31.9   122      2.89  84.3  90.1
str(study_t)
## tibble [12 × 2] (S3: tbl_df/tbl/data.frame)
##  $ month    : Factor w/ 12 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ mean_temp: num [1:12] -5.64 -5.78 -0.64 6.37 12.83 ...
ggplot() + 
geom_col(data = study_t, aes(x = month, y = mean_temp), fill = "red3", alpha = 0.7) + 
geom_line(data = study_p, aes(x = month, y = total_precip-80), col= "steelblue4", size = 1.5) +
  scale_y_continuous("Temperature (°C)", sec.axis = sec_axis(~ . + 80, name = "Precipitation (mm)")) +
    geom_ribbon(data = study_p, aes(x = month, ymin = (lower-80), ymax = (upper-80)), alpha = 0.3) +
  scale_x_discrete("Month") +
  theme_bw() +
   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

p.climograph <- ggplot() +
  geom_boxplot(data = pinery_temp_mean, aes(x = month, y = mean_temp_year, group = month), fill = "red3", alpha = 0.7) +
  geom_line(data = study_p, aes(x = month, y = total_precip-75), col= "steelblue4", size = 1) +
  geom_ribbon(data = study_p, aes(x = month, ymin = (lower-75), ymax = (upper-75)), fill = "steelblue4", alpha = 0.2) +
  scale_y_continuous("Temperature (°C)", sec.axis = sec_axis(~ . + 75, name = "Precipitation (mm)")) +
  scale_x_discrete("Month") +
  theme_bw() +
   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

p.climograph

#ggsave('p.climograph.pdf', p.climograph, units = 'cm', width = 20, height = 12)

5 Annual Data

london_temp_annual <- london_temp %>% filter(month == "Annual")
london_precip_annual <- london_precip %>% filter(month == "Annual")
pinery_temp_annual <- pinery_temp_mean %>% group_by(year) %>% summarise(mean_temp = mean(mean_temp_year))
## `summarise()` ungrouping output (override with `.groups` argument)
pinery_precip_annual <- pinery_precip_mean %>% group_by(year) %>% summarise(total_precip = sum(total_precip_year))
## `summarise()` ungrouping output (override with `.groups` argument)
london_temp_2 <- london_temp_annual %>% select(-month)
london_temp_2$location <- "London"
pinery_temp_annual$location <- "Pinery"
annual_temp <- rbind(london_temp_2,pinery_temp_annual) #join data sets

london_precip_2 <- london_precip_annual %>% select(-month)
london_precip_2$location <- "London"
pinery_precip_annual$location <- "Pinery"
annual_precip <- rbind(london_precip_2,pinery_precip_annual) #join data sets

Plot Temperature

annual_temp_plot <- ggplot() +
  # London Data
  geom_line(data=london_temp_annual, aes(x = year, y = mean_temp), colour = "lightseagreen", linetype = "dashed", alpha = 0.8) +
  # Pinery Data
  geom_line(data=pinery_temp_annual, aes(x = year, y = mean_temp), colour="lightsalmon4") +
  #Aesthetics
  xlab("Year") +
  ylab("Mean Annual Temperature (°C)") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  #ggtitle("Average Annual Temperature in \n Southwestern Ontario from 1883 to 2017") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

annual_temp_plot
## Warning: Removed 1 row(s) containing missing values (geom_path).

#ggsave('annual_temp.pdf', annual_temp_plot, units = 'cm', width = 15, height = 10)

Plot recent temperature

pinery_temp_annual_recent <- pinery_temp_annual %>% filter(year > 1949)

annual_temp_plot2 <- ggplot() +
  # Pinery Data
  geom_line(data=pinery_temp_annual_recent, aes(x = year, y = mean_temp)) +
  geom_point(data=pinery_temp_annual_recent, aes(x = year, y = mean_temp), size = 1)  +
  #Aesthetics
  xlab("Year") +
  ylab("Temperature (°C)") +
  scale_x_continuous(breaks=seq(1950,2020,10)) +
  #ggtitle("Average Annual Temperature in \n Southwestern Ontario from 1883 to 2017") + 
  scale_colour_manual(name="Location", breaks = c("London", "Pinery"),
                      values = c("blue","darkgreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
#ggsave('annual_temp.pdf', annual_temp_plot, units = 'cm', width = 20, height = 12)
annual_temp_plot2

Plot Precipitation

annual_precip_plot <- ggplot() +
  # London Data
  geom_line(data=london_precip_annual, aes(x = year, y = total_precip), colour = "lightseagreen", linetype = "dashed", alpha = 0.8) +
   # geom_smooth(data=london_precip_annual, aes(x = year, y = total_precip), se = TRUE, 
             # method = 'lm', size = 0.8, colour = "black", fill = "blue") +
  # Pinery Data
  geom_line(data=pinery_precip_annual, aes(x = year, y = total_precip, group = 1),  colour="lightsalmon4")+
   # geom_smooth(data=pinery_precip_annual, aes(x = year, y = total_precip), se = TRUE, 
            #  method = 'lm', size = 0.8, colour = "black", fill = "darkgreen") +
  #Aesthetics
  xlab("Year") +
  ylab("Total Annual Precipitation (mm)") +
  scale_x_continuous(breaks=seq(1880,2020,20)) +
  #ggtitle("Annual Precipitation in Southwestern Ontario \nfrom 1883 to 2017") + 
  scale_colour_manual(name="Location", breaks = c("London", "Pinery"),
                      values = c("blue","darkgreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

annual_precip_plot
## Warning: Removed 17 row(s) containing missing values (geom_path).

#ggsave('annual_precip.pdf', annual_precip_plot, units = 'cm', width = 15, height = 10)
annual_plot <- grid.arrange(annual_temp_plot, annual_precip_plot, ncol= 2)
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 17 row(s) containing missing values (geom_path).

#ggsave('annual_plot.pdf', annual_plot, units = 'cm', width = 23, height = 8)

6 Seasonal Data

winter_t <- pinery_temp %>% filter(month == 1 | month == 2 | month == 3)
winter_t$season <- "winter"
spring_t <- pinery_temp %>% filter(month == 4 | month == 5 | month == 6)
spring_t$season <- "spring"
summer_t <- pinery_temp %>% filter(month == 7 | month == 8 | month == 9)
summer_t$season <- "summer"
fall_t <- pinery_temp %>% filter(month == 10 | month == 11 | month == 12)
fall_t$season <- "fall"
pinery_temp <- rbind(winter_t, spring_t, summer_t, fall_t)
pinery_temp$season <- as.factor(pinery_temp$season) #Make a factor 
str(pinery_temp)
## 'data.frame':    1944 obs. of  5 variables:
##  $ year     : int  1885 1885 1885 1886 1886 1886 1887 1887 1887 1888 ...
##  $ month    : int  1 2 3 1 2 3 1 2 3 1 ...
##  $ mean_temp: num  -9.4 -13.5 -8.9 -7.2 -6.9 -1 -8.4 -5.4 -3.8 -9.6 ...
##  $ station  : Factor w/ 9 levels "Arkona","Centralia A",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ season   : Factor w/ 4 levels "fall","spring",..: 4 4 4 4 4 4 4 4 4 4 ...
pinery_temp_station <- pinery_temp %>% 
 group_by(year, season, station) %>% 
 summarise(mean_temp = mean(mean_temp))
## `summarise()` regrouping output by 'year', 'season' (override with `.groups` argument)
pinery_temp_season <- pinery_temp_station %>%
  group_by(year, season) %>%
  summarise(mean_temp_season = mean(mean_temp) %>% round(2), sd_temp = sd(mean_temp) %>% round(2), n = n()) %>% 
  mutate(se_temp = (sd_temp / sqrt(n)) %>% round(2), lower = mean_temp_season - se_temp, upper = mean_temp_season + se_temp) 
## `summarise()` regrouping output by 'year' (override with `.groups` argument)
head(pinery_temp_season)
## # A tibble: 6 x 8
## # Groups:   year [2]
##    year season mean_temp_season sd_temp     n se_temp lower upper
##   <int> <fct>             <dbl>   <dbl> <int>   <dbl> <dbl> <dbl>
## 1  1885 fall               2.7       NA     1      NA    NA    NA
## 2  1885 spring            10.9       NA     1      NA    NA    NA
## 3  1885 summer            17.6       NA     1      NA    NA    NA
## 4  1885 winter           -10.6       NA     1      NA    NA    NA
## 5  1886 fall               1.33      NA     1      NA    NA    NA
## 6  1886 spring            12.4       NA     1      NA    NA    NA
winter_p <- pinery_precip %>% filter(month == 1 | month == 2 | month == 3)
winter_p$season <- "winter"
spring_p <- pinery_precip %>% filter(month == 4 | month == 5 | month == 6)
spring_p$season <- "spring"
summer_p <- pinery_precip %>% filter(month == 7 | month == 8 | month == 9)
summer_p$season <- "summer"
fall_p <- pinery_precip %>% filter(month == 10 | month == 11 | month == 12)
fall_p$season <- "fall"
pinery_precip <- rbind(winter_p, spring_p, summer_p, fall_p)
pinery_precip$season <- as.factor(pinery_precip$season) #Make a factor 

pinery_precip_station <- pinery_precip %>%
  group_by(year, season, station) %>%
  summarise(total_precip = mean(total_precip))
## `summarise()` regrouping output by 'year', 'season' (override with `.groups` argument)
pinery_precip_season <- pinery_precip_station %>% 
  group_by(year, season) %>% 
  summarise(total_precip_season = mean(total_precip) %>% round(2), sd_precip = sd(total_precip) %>% round(2), n = n()) %>%
  mutate(se_precip = sd_precip / sqrt(n), lower = total_precip_season - sd_precip, upper = total_precip_season + sd_precip %>% round(2)) 
## `summarise()` regrouping output by 'year' (override with `.groups` argument)
head(pinery_precip_season)
## # A tibble: 6 x 8
## # Groups:   year [2]
##    year season total_precip_season sd_precip     n se_precip lower upper
##   <int> <fct>                <dbl>     <dbl> <int>     <dbl> <dbl> <dbl>
## 1  1883 fall                  69.3     NA        1    NA      NA    NA  
## 2  1883 spring               101.      NA        1    NA      NA    NA  
## 3  1883 summer                83.8     NA        1    NA      NA    NA  
## 4  1883 winter                65.2     NA        1    NA      NA    NA  
## 5  1884 fall                  91.5      1.01     2     0.714  90.5  92.5
## 6  1884 spring                49.1      1.51     2     1.07   47.6  50.6
winter_tl <- london_temp %>% filter(month == "Jan" | month == "Feb" | month == "Mar")
winter_tl$season <- "winter"
spring_tl <- london_temp %>% filter(month == "Apr" | month == "May" | month == "Jun")
spring_tl$season <- "spring"
summer_tl <- london_temp %>% filter(month == "Jul" | month == "Aug" | month == "Sep")
summer_tl$season <- "summer"
fall_tl <- london_temp %>% filter(month == "Oct" | month == "Nov" | month == "Dec")
fall_tl$season <- "fall"
london_temp <- rbind(winter_tl, spring_tl, summer_tl, fall_tl)
london_temp$season <- as.factor(london_temp$season) #Make a factor 
str(london_temp)
## 'data.frame':    1620 obs. of  4 variables:
##  $ year     : int  1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 ...
##  $ month    : Factor w/ 17 levels "Annual","Apr",..: 7 7 7 7 7 7 7 7 7 7 ...
##  $ mean_temp: num  NA -10.1 -9.2 -7.2 -8.3 -10.1 -3 -0.7 -4.3 -7.9 ...
##  $ season   : Factor w/ 4 levels "fall","spring",..: 4 4 4 4 4 4 4 4 4 4 ...
london_temp_season <- london_temp %>% 
  group_by(year, season) %>% 
  summarise(mean_temp_season = mean(mean_temp))
## `summarise()` regrouping output by 'year' (override with `.groups` argument)
winter_pl <- london_precip %>% filter(month == "Jan" | month == "Feb" | month == "Mar")
winter_pl$season <- "winter"
spring_pl <- london_precip %>% filter(month == "Apr" | month == "May" | month == "Jun")
spring_pl$season <- "spring"
summer_pl <- london_precip %>% filter(month == "Jul" | month == "Aug" | month == "Sep")
summer_pl$season <- "summer"
fall_pl <- london_precip %>% filter(month == "Oct" | month == "Nov" | month == "Dec")
fall_pl$season <- "fall"
london_precip <- rbind(winter_pl, spring_pl, summer_pl, fall_pl)
london_precip$season <- as.factor(london_precip$season) #Make a factor 
str(london_precip)
## 'data.frame':    1620 obs. of  4 variables:
##  $ year        : int  1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 ...
##  $ month       : Factor w/ 17 levels "Annual","Apr",..: 7 7 7 7 7 7 7 7 7 7 ...
##  $ total_precip: num  NA 106.2 96.1 142.5 86.9 ...
##  $ season      : Factor w/ 4 levels "fall","spring",..: 4 4 4 4 4 4 4 4 4 4 ...
london_precip_season <- london_precip %>% 
  group_by(year, season) %>% 
  summarise(total_precip_season = mean(total_precip))
## `summarise()` regrouping output by 'year' (override with `.groups` argument)

Make a factor

london_temp_season$season <- factor(london_temp_season$season, levels = c("winter", "spring", "summer", "fall"))
pinery_temp_season$season <- factor(pinery_temp_season$season, levels = c("winter", "spring", "summer", "fall"))

6.1 Combine London & Pinery Records

pinery_temp_season$location <- "Pinery"
london_temp_season$location <- "London"
season_temp <- rbind(pinery_temp_season, london_temp_season)

pinery_precip_season$location <- "Pinery"
london_precip_season$location <- "London"
season_precip <- rbind(pinery_precip_season, london_precip_season)

season_temp$location <- factor(season_temp$location, levels = c("Pinery", "London"))
season_precip$location <- factor(season_precip$location, levels = c("Pinery", "London"))

Plot Temperature

w_temp <- ggplot() +
  geom_line(data=(season_temp %>% filter(season == "winter")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature (°C)") + ggtitle("Winter") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

sp_temp <- ggplot() +
  geom_line(data=(season_temp %>% filter(season == "spring")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature (°C)") + ggtitle("Spring") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
    scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

su_temp <- ggplot() +
  geom_line(data=(season_temp %>% filter(season == "summer")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature (°C)") + ggtitle("Summer") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
    scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

fa_temp <- ggplot() +
  geom_line(data=(season_temp %>% filter(season == "fall")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature (°C)") + ggtitle("Fall") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
    scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),  legend.position = "none") 

seasonal_temp_plot <- grid.arrange(w_temp, sp_temp, su_temp, fa_temp, nrow = 2, ncol = 2)
## Warning: Removed 1 row(s) containing missing values (geom_path).

#ggsave('seasonal_temperature.pdf', seasonal_temp_plot, units = 'cm', width = 20, height = 12)

Plot Precipitation

w_precip <- ggplot() +
  geom_line(data=(season_precip %>% filter(season == "winter")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Winter") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
      scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

sp_precip <- ggplot() +
  geom_line(data=(season_precip %>% filter(season == "spring")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Spring") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
      scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

su_precip <- ggplot() +
  geom_line(data=(season_precip %>% filter(season == "summer")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Summer") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
      scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

fa_precip <- ggplot() +
  geom_line(data=(season_precip %>% filter(season == "fall")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Fall") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
      scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),  legend.position = "none") 

seasonal_precip_plot <- grid.arrange(w_precip, sp_precip, su_precip, fa_precip, nrow = 2, ncol = 2)
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 15 row(s) containing missing values (geom_path).

## Warning: Removed 15 row(s) containing missing values (geom_path).
## Warning: Removed 16 row(s) containing missing values (geom_path).

#ggsave('seasonal_precipitation.pdf', seasonal_precip_plot, units = 'cm', width = 20, height = 12)

6.2 Temperature Anomaly

https://tamino.wordpress.com/2018/07/29/why-use-temperature-anomaly/

Average defined as 1951-1980 (NASA standard)

pinery_average <- season_temp %>% filter(location == "Pinery") %>% filter(year > 1950 & year < 1981)
pinery_average <- pinery_average %>% group_by(season) %>% summarise(pinery_average = mean(mean_temp_season))
## `summarise()` ungrouping output (override with `.groups` argument)
london_average <- season_temp %>% filter(location == "London") %>% filter(year > 1950 & year < 1981)
london_average <- london_average %>% group_by(season) %>% summarise(london_average = mean(mean_temp_season))
## `summarise()` ungrouping output (override with `.groups` argument)
pinery_anomaly <- right_join(season_temp, pinery_average)
## Joining, by = "season"
pinery_anomaly <- mutate(pinery_anomaly, anomaly = mean_temp_season - pinery_average) %>% filter(location == "Pinery")

london_anomaly <- right_join(season_temp, london_average)
## Joining, by = "season"
london_anomaly <- mutate(london_anomaly, anomaly = mean_temp_season - london_average) %>% filter(location == "London")

season_anomaly <- full_join(pinery_anomaly, london_anomaly) %>% select(-pinery_average, -london_average)
## Joining, by = c("year", "season", "mean_temp_season", "sd_temp", "n", "se_temp", "lower", "upper", "location", "anomaly")
w_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "winter")), aes(x = year, y = anomaly, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Winter") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

sp_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "spring")), aes(x = year, y = anomaly, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Spring") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

su_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "summer")), aes(x = year, y = anomaly, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Summer") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

fa_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "fall")), aes(x = year, y = anomaly, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Fall") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),  legend.position = "none") 

seasonal_anomaly_plot <- grid.arrange(w_anomaly, sp_anomaly, su_anomaly, fa_anomaly, nrow = 2, ncol = 2)
## Warning: Removed 1 row(s) containing missing values (geom_path).

#ggsave('seasonal_anomaly.pdf', seasonal_anomaly_plot, units = 'cm', width = 20, height = 12)

Plot recent temperature anomaly

w_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "winter", year > 1949, location == "Pinery")), aes(x = year, y = anomaly)) +
  geom_point(data=(season_anomaly %>% filter(season == "winter", year > 1949, location == "Pinery")), aes(x = year, y = anomaly), size = 0.5) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Winter") +
  scale_x_continuous(limits = c(1950, 2020), breaks=seq(1950,2020,10)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

sp_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "spring", year > 1949, location == "Pinery")), aes(x = year, y = anomaly)) +
 geom_point(data=(season_anomaly %>% filter(season == "spring", year > 1949, location == "Pinery")), aes(x = year, y = anomaly), size = 0.5) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Spring") +
  scale_x_continuous(limits = c(1950, 2020), breaks=seq(1950,2020,10)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

su_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "summer", year > 1949, location == "Pinery")), aes(x = year, y = anomaly)) +
  geom_point(data=(season_anomaly %>% filter(season == "summer", year > 1949, location == "Pinery")), aes(x = year, y = anomaly), size = 0.5) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Summer") +
  scale_x_continuous(limits = c(1950, 2020), breaks=seq(1950,2020,10)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

fa_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "fall", year > 1949, location == "Pinery")), aes(x = year, y = anomaly)) +
    geom_point(data=(season_anomaly %>% filter(season == "fall", year > 1949, location == "Pinery")), aes(x = year, y = anomaly), size = 0.5) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Fall") +
  scale_x_continuous(limits = c(1950, 2020), breaks=seq(1950,2020,10)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),  legend.position = "none") 

seasonal_anomaly_plot <- grid.arrange(w_anomaly, sp_anomaly, su_anomaly, fa_anomaly, nrow = 2, ncol = 2)

7 Save Outputs

#write.csv(pinery_temp_mean, "outputs/monthly_temperature.csv")
#write.csv(pinery_precip_mean, "outputs/monthly_precipitation.csv")
#write.csv(season_temp, file = "outputs/season_temperature.csv")
#write.csv(season_precip, file = "outputs/season_precipitation.csv")
#write.csv(pinery_temp_annual, file = "outputs/annual_temp.csv")

8 Reproducibility

Sys.time()
## [1] "2020-08-05 16:30:07 PDT"
git2r::repository()
## Local:    master /Users/JenBaron/Documents/UWO/UWO NSERC/Statistical Analysis/Climate/Climate_Pinery
## Remote:   master @ origin (https://github.com/JenBaron/Climate_Pinery.git)
## Head:     [a85d3a9] 2020-08-05: Fill missing values with London & ClimateNA
sessionInfo()
## R version 4.0.0 (2020-04-24)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.5
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gridExtra_2.3   GGally_1.5.0    forcats_0.5.0   stringr_1.4.0  
##  [5] dplyr_1.0.0     purrr_0.3.4     readr_1.3.1     tidyr_1.0.2    
##  [9] tibble_3.0.1    tidyverse_1.3.0 ggplot2_3.3.1  
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.0   xfun_0.14          haven_2.2.0        lattice_0.20-41   
##  [5] colorspace_1.4-1   vctrs_0.3.1        generics_0.0.2     htmltools_0.4.0   
##  [9] yaml_2.2.1         utf8_1.1.4         rlang_0.4.6        pillar_1.4.4      
## [13] glue_1.4.1         withr_2.2.0        DBI_1.1.0          RColorBrewer_1.1-2
## [17] dbplyr_1.4.3       modelr_0.1.6       readxl_1.3.1       plyr_1.8.6        
## [21] lifecycle_0.2.0    munsell_0.5.0      gtable_0.3.0       cellranger_1.1.0  
## [25] rvest_0.3.5        evaluate_0.14      labeling_0.3       knitr_1.28        
## [29] fansi_0.4.1        broom_0.5.6        Rcpp_1.0.4.6       scales_1.1.1      
## [33] backports_1.1.7    jsonlite_1.6.1     farver_2.0.3       fs_1.4.1          
## [37] hms_0.5.3          digest_0.6.25      stringi_1.4.6      grid_4.0.0        
## [41] cli_2.0.2          tools_4.0.0        magrittr_1.5       crayon_1.3.4      
## [45] pkgconfig_2.0.3    ellipsis_0.3.1     xml2_1.3.2         reprex_0.3.0      
## [49] lubridate_1.7.8    reshape_0.8.8      assertthat_0.2.1   rmarkdown_2.1     
## [53] httr_1.4.1         rstudioapi_0.11    R6_2.4.1           git2r_0.26.1      
## [57] nlme_3.1-147       compiler_4.0.0